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Abstract — This paper presents a trajectory generation 
method that optimizes a quadratic cost functional with respect 
to linear system dynamics and to linear input and state con- 
straints. The method is based on continuous-time flatness-based 
trajectory generation, and the outputs are parameterized using 
a polynomial basis. A method to parameterize the constraints 
is introduced using a result on polynomial nonpositivity. The 
resulting parameterized problem remains linear-quadratic and 
can be solved using quadratic programming. The problem can 
be further simplified to a linear programming problem by 
linearization around the unconstrained optimum. The method 
promises to be computationally efficient for constrained systems 
with a high optimization horizon. As application, a predictive 
torque controller for a permanent magnet synchronous motor 
which is based on real-time optimization is presented. 

I. Introduction 

Trajectory optimization in real-time is an important part 
of many modern control systems, for instance, predictive 
control. The trajectory optimization problem must be solved 
in the sampling interval, determined by the plant dynamics. 
A computationally efficient solution, however, encounters 
two major obstacles: first, the optimization horizon has a 
minimum length, either to avoid suboptimality or as stability 
criterion, and second, the trajectory must satisfy input and 
state constraints to be feasible. For fast systems, many of 
the existing constrained predictive control schemes, which 
are based on numerical iterations, are not applicable. 

Concerning the horizon length problem, the continuous 
approach to predictive control is of interest [1]. Using differ- 
ential flatness, the optimal control problem can be rewritten 
as output optimization problem. The basis function approach 
is applied to obtain a finite-parameter optimization problem 
[2] [3] [4] [5], analogeously to discrete-time optimization. 
A long optimization horizon is, however, obtained with 
comparably few optimization parameters. 

A remaining issue regarding computational efficiency is 
the inclusion of constraints. The classical method is the 
application of penalty functions [3]. As an alternative, a 
coordinate transformation was proposed to modify the con- 
strained problem to a nonlinear unconstrained problem [6], 
solvable with nonlinear calculus of variations methods. The 
existing optimization methods with penalty functions or 
nonlinear coordinate transformations lead to a nonlinear or 
non-linear-quadratic problem, increasing the computational 
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burden. If only feasibility is of interest, numerical procedures 
applying time scaling to slow down some variables can be 
used [1]. 

This paper treats the special case of the linearly con- 
strained linear-quadratic problem. The linear structure of the 
system is exploited in the transformation of the problem to a 
finite-parameter optimization problem. The linear-quadratic 
structure of the problem is maintained. The unconstrained 
problem is solved algebraically as it is convex. A quadratic 
or linear programming (QP/LP) solver is then applied for 
constraint handling to modify the solution to a feasible 
trajectory. The solution is computationally efficient, as first 
a continuous parameterization is performed, and second a 
linear programming solver is used for the quadratic problem 
by linearizing around the unconstrained optimum. To the 
knowledge of the authors, QP and LP have so far not been 
applied to this type of continuous problems. 

The paper is organized as follows. Section II defines the 
optimization problem and presents inversion-based linear- 
quadratic trajectory optimization using basis functions, with 
the special case of power series. Section III proposes the 
novel constraint handling method, first using a result on 
polynomial nonnegativity, second rewriting the problem as 
linear programming problem. Section IV shows numerical 
results of a synchronous motor predictive torque controller, 
which is a multivariable system with a high sampling rate and 
both input and state constraints. The algorithm is computable 
in real time: 2 ms prediction are implemented at a sampling 
rate of 10 kHz. The experimental results are presented in 
[21]. 

ii. preliminaries: continuous flatness-based 
Linear-Quadratic Trajectory Optimization 

A. Optimal Control Problem 

The optimal control problem considered is to generate tra- 
jectories on a finite time horizon T minimizing the quadratic 
cost functional 



(x T (t)Qx(t) +u T (i)Ru(t)) dt 
+ (x(T)-x*) T P(x(T)-x*) 



(1) 



with states x(t) E M™, inputs u(t) e 



p e R nx ", Q g 

final state x* e K 



weight matrices 
nxn and R e R mxm , and the desired 
The weight matrices are assumed to be 
positive definite and symmetric. Any quadratic functional is 
eligible, for instance the reference x* can also be included 
in the cost integral. Quadratic cost functionals are preferred 



in predictive control as they will provide good closed-loop 
behavior [7]. They can also describe physical costs better 
than other norms. Furthermore, if the weight matrices are 
positive definite, the optimization problem is convex and 
easily solvable: a unique optimum exists and is found by 
solving first-order necessary conditions. The states and inputs 
are constrained to the time invariant linear dynamics of the 
multi-input multi-output system 

x(t) = Ax(t) + Bu(t) (2) 

with system matrix A G R nxn and input matrix B G 
jgmxm S y S t em [ s further assumed to be controllable. The 
optimization is subject to N c linear input and state constraints 

G x x(t) + G u u(t) + g < 0, V* G [0, T] (3) 

with G x G R N °* n , G u G R N °* m , g G R N ", as well as the 
n initial conditions 

x(0) = x . (4) 

Terminal constraints could be imposed too, but may yield an 
unfeasible problem in the presence of input or state bounds. 

B. Parameterization of the System Variables using Flatness 

By definition, a linear system is said to be differentially 
flat if there exist m output functions 

y f (t) = C f x(t), (5) 

with y f (t) G R"\ Cf G R mx ", such that all states x(t) 
and inputs u(i) can be expressed as linear combination of 
the flat outputs and a finite number of their derivatives [8] 
[9]. Differential flatness can be interpreted as transformation 
into controller canonical form, and in the linear case, it is 
equivalent to controllability [9]. The flat outputs yy(t) are 
thus the controller canonical form outputs, and the canonical 
form state vector is 

Z <X> = ((y/i>y/i>">2//i 1-1) ) , (yf2,yj2,-,y% 2 ' 1) ) , 

•-, (y/m,y/m,-,y/ r ^ _1) )) (6) 

with z(t) G R™ and where 7"$ is the corresponding vector 
relative degree with regard to the output yji (yfi being the 
i-th flat output resp. z-th element of in (5), see [9]). Note 
that 53i=i r i = n - The differential parameterization of the 
system variables is [8] [9] 

* = s *(y/ J y/.".y?" 1) ) s ( ? ) 

u = S u (y / ,y / ,.,yf) ) (8) 

where k E R is k = max{ri}. The functions S x and a u 
are linear functions of z resp. of z and z, as they are derived 
from the controller canonical form transformation. See that 
(O is the inverse transformation from controller canonical 
form, and (0 follows from (0 and (0. 



C. Output Parameterization with a Polynomial Basis 

The Ritz-Galerkin method, also called basis function ap- 
proach, is a direct method to find an approximate solution 
to an optimal control problem [10]. In the flatness-based ap- 
proach, the outputs are parameterized as a linear combination 
of time-variant basis functions. It is comparable to the control 
parameterization Ritz method (CPRM) [11], where the inputs 
u are parameterized. Convergence of this method (combined 
with linear splines) was studied in [12] and more generally 
in [13]. 

As basis functions, polynomials are chosen. Their simplic- 
ity is exploited for the further developments, main advantage 
being that the parameters enter linearly, for instance they 
allow to prove that the transformed cost remains convex. 
They also allow constraint transformation based on a result 
on polynomial nonpositivity in section III. Furthermore, 
the dynamic shape of the resulting trajectory can be well 
analyzed. 

Higher-order polynomials (for instance Laguerre and Leg- 
endre) have been used in many applications, but here, 
power series, the simplest form of polynomials, are used. 
The methods and results are all applicable to higher-order 
polynomials, there, the undetermined coefficients still enter 
linearly and they might inherit numerical advantages (power 
series are numerically stable only up to degree 12.. 15), but 
the notation would be less comprehensible. In [14] it was 
shown that only the polynomial degree but not the type of 
series are important for convergence. The choice is a simple 
approximate, but yet, it turns out to be quite accurate. 

The system outputs are defined as degree N power series 

N P 

= Y, a « jj ' l = 1 - m > 1 G [°> T l (9) 

with aij G R. This is not an approximation of a system step 
response, but merely an assumption that the optimal solution 
can be described as a polynomial. The system response is 
exact, only the output is reduced in dimensionality. 

The original problem is thereby transformed to a finite- 
parameter optimization problem where a set of constants is 
searched. Fig. Q] demonstrates the computational advantage. 
For instance, in the example in section IV, the prediction 
horizon is 2 ms (the slowest time constant in the system) 
and the sampling rate is 10 kHz (determined by the power 
stage). A discrete description takes 20 parameters, but here, 
only 6 are sufficient to describe a setpoint change. A more 
complicated trajectory may not be expected on most plants. 

The n initial conditions define n of the parameters ay, 
they are computed via the transformed state z(0) in © (See 
that the initial time is to = 0, such that only the first element 
of the power series is nonzero, same for the derivatives). For 
further notational simplicity, the vector a is defined as the 
vector of the undetermined coefficients 

« = (airi! ■;OLljf,a2r 2 , -,OL2N, - j OW m , •■, CK m N) T , (10) 

with a £ W N with N' — m x N — n, and the indices 
n represent the respective vector relative degree of the flat 
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Fig. 1. Trajectories with 6 free parameters. Discrete-time horizon is 
maximum 0.6 ms, continuous trajectory is well-conditioned at the desired 
horizon length of 2 ms. 

output y fi . 

Based on the flatness parameterization © and (0, and 
knowing that t are polynomials, the system variables 
can be parameterized with the coefficients a. yielding the 
functions 

x(t) = r x (cM), (ii) 

u(t) = r tt (a,i). (12) 

These functions are degree N polynomials of time t, where 
the coefficients a. enter linearly. 

It is seen that N > maxjr,;}, i = l..m, is mandatory 
to have degrees of freedom on the trajectory. It should be 
remarked that when choosing a higher N, the error bounds 
of the basis function approach become smaller and the 
dynamical response increases, at the cost of higher com- 
putational requirements and the risk of less good numerical 
conditioning of the parameters. 

D. Remark: Generality of the Parameterization for Linear 
Systems with Power Series Outputs 

It is remarked that the same parameterization (fTTT i and (fT2l 
can be obtained for arbitrary non-flat outputs of a controllable 
and observable linear system, if the output is defined as 
polynomial series with undetermined coefficients [15]. Then, 
the outputs do not need to be redefined to flat, or controller 
canonical form outputs. 

E. Parameterizing the Cost Functional 

The parameterization of the system variables (fTTT i. (fT2l 
is applied to transform the cost functional (UJi to a finite- 
parameter cost function. 

Proposition II. 1 

Conditioning the cost functional. 

The cost functional (T7J with the substitutions ( 1771 ) and ( 1721 ) 
yields a biaffine function in the undetermined parameters a. 
of the type 

J (a) = a T Ka+k T a + fc (13) 

with k £ R, k £ R N ' and K £ R N ' xN '. | 

The proof follows from a straightforward computation 
with the knowledge that ( fTTT i. (fT2l are polynomials with 
affine coefficients a, and is omitted. 



Eq. dT3l i is derived directly from (1) by substitution of 
(7), (8) and (9) (or (11) and (12)). Practical (symbolical) 
calculations are performed using a computer algebra tool, 
matrix K is the Hessian matrix of J, and vector k follows 
from the gradient of J. 

Proposition II.2 

Convexity of the conditioned cost function. 
The conditioned cost function ( 1731 ) is convex over the set 
a £ M. N if matrices Q, R and P in (fJ]) are positive definite 
and real-symmetric matrices. | 

This statement of convexity is not obvious, as the trans- 
formed cost function ( fT3l has a higher dimension as the orig- 
inal cost functional ([T]) (TV 7 instead of n). It is proven for the 
unconstrained linear-quadratic case. This proof follows from 
the knowledge that the states and inputs are polynomials 
with coefficients linearly dependent on a, and is sketched 
in appendix I. As only linear constraints are regarded, the 
result can be assumed valid also for the constrained case. 

This result allows to apply first-order conditions of op- 
timality, and will guarantee a solution in finite time as an 
unique minimum of J exists. 

III. Main Results: Constraint Handling 

The previous section presented the unconstrained opti- 
mization of a controllable linear system with a flat output pa- 
rameterized with a polynomial basis. This section introduces 
a computationally efficient method for including constraints 
in trajectory generation. It is a follow-up to the previous 
section. 

A. Parameterizing the Constraints 

The input and state constraints (O are transformed with 
the system variable parameterizations (fTTl) and (fT2"l) . With 
the knowledge that the parameterizations ( fTTl ) and ( fT2l ) are 
univariate polynomials in t, each of the constraints in (O is 
also a polynomial in t of order N where the coefficients are 
affine functions in a. Thus, the constraints are rewritten as 

N 

Pk(t) = (9kio +£ki «) < 0- k = l - N ct £ [0,T], 

(14) 

with gfcjQ £ K and g ki £ M w . Checking constraints 
over a time interval t £ [0,T] is difficult. The typical 
solutions for such continuous problems are penalty functions 
in the cost functional, or nonlinear state transformations. 
In the following, the constraints are transformed to allow 
a direct check on a. independent of t, while maintaining the 
simplicity of a linear constraint. 

Proposition III.l 

Sufficient affine conditions for the constraints. 

The univariate polynomials Pk (t) are positive on the finite 



interval t G [0, T] if the conditions 

P fe (0)<0, Vk = l..N c , (15) 




-A-P fe (0)<0, Vp = l..N,k = l..N c , (16) 



with a constant A G M > are satisfied. | 

The constraints, written in the introduced notations, are 
(5fcoo+gM ") <0, fc = l.JV c (17) 

and 

N / rpi \ 

( 9m + gLa)—p l - A ■ ( 5fc00 + gl a J < 0, 

i=0 ^ ' 

p = l..N,k= l.JV c . (18) 

The underlying idea is shown in Fig. [2] If a polynomial 
trajectory P(t) of degree N is constrained at A^ + l sampling 
points on the interval [0, T] to be nonpositive, then, in the 
worst case, this polynomial reaches a maximum of P(t) = 
—A • .P(O). This maximum offset of —A • P(0) is applied 
as interleaf to the constraint boundary (the second term in 
(TToTl). The proof of this proposition is presented in appendix 
II, where also the value of the constant A is found, which is a 
constant depending on the polynomial degree N only. After 
transformation of the constraints, they can be used directly 
in an optimization procedure, as the result is an affine and 
purely parametric constraint. They are affine functions of a. 
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Fig. 2. A degree N polynomial trajectory (black) constrained at TV + 1 
points to be nonpositive (red arrows) will not exceed the upper bound of 
—A ■ P(0) (green line). Here, N = 6, P(0) = -1 and A = 0.037. 

The resulting set of constraints is sufficient, but not neces- 
sary, thus they guarantee maintainance of the constraints, but 
are too restrictive. The offset between the necessary and the 
sufficient conditions is the factor A, and the restriction is, for 
example, 6.4% for N = 3, 1.2% for N = 10, and even less 
for higher N. It is assumed that this restriction is acceptable 
for most applications regarding the inherited computational 
advantages. 

Exact results on necessary and sufficient conditions on 
positivity of univariate polynomials over an interval [0, T] 
exist. Semidefinite programming (SDP) is applied to estab- 
lish linear matrix inequalities (LMI), which, just as for the 
presented results, are independent of t and can directly be 



used in a program [16]. These methods, however, require 
a numerical iterative algorithm to establish the LMIs, and 
can therefore not be considered if computational efficiency is 
important. Due to updated initial conditions, the constraints 
must be reconditioned in each trajectory generation cycle, 
and the computational burden adds up. 

B. Trajectory Generation in Finite Time 

In the previous section, the optimal control problem has 
been transformed to a finite-parameter optimization problem. 
The transformed problem can be solved in finite time using 
quadratic programming (QP) to find the parameters a, as 
the cost function is quadratic in a and as the constraints 
are affine functions of a. This is an interesting result for 
continuous systems, as QP is a standard method in discrete- 
time optimization. Furthermore, the number of parameters 
is decoupled from the horizon length, thus higher prediction 
horizons can be reached with less parameters. 

An even more convincing advantage of the presented 
method shall be a further reduction of the computational 
complexity. In the following, the problem will be approx- 
imated such that it is solvable using linear programming 
(LP) techniques. The approximation yields a feasible solu- 
tion (thus it is satisfying the constraints) with a bounded 
suboptimality. 

C. Transformation to a Least Distance Problem 

In the first step, a transformation is performed to obtain a 
simpler cost functional. 

Proposition III.2 

Reformulation as least-distance problem. 
The cost functional can be transformed to 

J=ff+c, (19) 

where f — (/i, /2, •■, /w) T € is the vector of trans- 

formed parameters and c G K a constant. | 

The new parameter vector / is defined as 

/ = F(a - am), (20) 

where ap G M. N is the unconstrained minimum of the 
parameterized cost functional (fT3l ) 

ao = -^K _1 k, (21) 

and matrix F G M. N xN is defined such that 

F T F = K (22) 

with K from eq. dT3T >. Computation is not an issue as K is 
positive definite and symmetric, for instance, in the notations 
in appendix I one would set F = BA. If this decomposition 
is not known, F can be computed via the Cholesky matrix 
decomposition F T = cholesky (K T ). 

Transformation < f20b must also be applied to transform 
the constraints. This is done by substituting the inverse 



transformation 

a = a Q + F- 1 / (23) 

to the constraint equation, yielding directly the affine con- 
straints in terms of /. The same equation is applied to 
retransform the obtained results into the original parame- 
ters. Contrary to the amount of parameters, the amount of 
constraints is not increased. 

The least-distance problem is simpler to solve as the orig- 
inal QP problem, as some computations become obsolete in 
the iterations [17]. The transformation renders the quadratic 
programming more efficient and increases numerical stabil- 
ity. 

D. Transformation to a Linear Programming Problem 

Now the least distance problem can be approximated for a 
solution using linear programming. The L 2 -norm is rewritten 
as Li-norm. The variables in linear programs are limited to 
positive numbers, thus the variables / are replaced by 

fi — fip fin ) (24) 

with f ip , f in e M and f ip > 0, f in > 0. 
Proposition III.3 

Linear Cost Function for the Least-Distance Problem. 
The cost function approximation 

N N N 

i—1 i—1 i—1 

(25) 

with fip, fi n > 0, yields a feasible solution with bounded 
suboptimality for a least-distance problem. | 

This approximation implies a large offset between the 
linearized cost J' and the correct cost J. However, this is not 
of importance, as the goal is to find a point in the parameter 
space / that is feasible and least-distance to the origin. Under 
this aspect, the suboptimality is not the variation of the cost, 
but of the difference of the found parameters in the squared 
distance to the origin. The contour lines of the quadratic cost 
are circles around the origin, whereas these of the linear cost 
are lozenges [7]. 

It should be remarked that |/j | is not necessarily equivalent 
to fip + fi n , but in the reverse, if fa p + fi n is minimized, it 
is equivalent to \fi \ as f p , fi n > 0, such that at least one of 
each f^, f in is zero. 

The advantages of a linear programming solver are ob- 
vious, increased reliability and reduced computational com- 
plexity, and the computational burden only grows linearly 
with the degree of the output polynomial N. It represents a 
significant computational saving compared to when applying 
QP as solver. 

If no constraint is active, the solution is exact. If a con- 
straint is active, the worst case is when the active constraint 
vertex is parallel to a contour line (i.e. parallel to a side of 



the lozenge). It can be shown that the resulting cost is 

J' = J a + N' x J c (26) 

in the worst case, where J' is the cost with the linear 
program, Jo the cost of the unconstrained problem, Jc 
the extra cost when considering constraints in a quadratic 
program, and TV' = dim(a) the amount of free parameters 
a. Thus, the linear programming solver yields the worst- 
case suboptimality J' /(Jo + Jc) of up to N' times the 
cost. If not all constraints are active, or if the constraints 
intersect the contour lines, suboptimality will be considerably 
less. The bound of the suboptimality and the consideration 
that the overall suboptimality is relative to the cost of the 
unconstrained problem makes the simplification of a linear 
programming solver acceptable for many applications. 

IV. Example: Predictive Torque Control of 
Permanent-Magnet Synchronous Machines 

To present the advantages and the good functionality of 
the presented trajectory optimization scheme, a predictive 
torque controller for a permanent-magnet synchronous ma- 
chine (PMSM) is presented. So far, this plant has been 
controlled via generalized predictive control (GPC), via finite 
control set MPC (FS-MPC) and via the explicit solution [18]. 
Even though the methods provided good results, they were 
either unconstrained, or limited to few prediction steps. The 
presented method obtains a high optimization horizon (2 ms) 
while respecting all constraints. 

Its suboptimality in 1) the assumption of a polynomial 
output, 2) the restrictive constraint parameterization (of 2%), 
and 3) the use of LP instead of QP is analyzed. 

A cascaded control structure is chosen [19]. The speed 
loop is controlled via a PI controller, and the current (and 
torque) loop is the predictive controller. The electrical sub- 
system of a non-salient PMSM, consisting of the torque- 
generating and field-generating currents i q and id (peak 
values), and the machine torque tm as output, is given as 

d R 1 

-^id = -jrid + n p LJMiq + jrV d (27) 

d . R . „ 1 

— i q = -j^h - n p tu M id - n p uj M K + —v q (28) 

3 

tm = ^n p Ki q (29) 

which is linearized by the assumption that speed is constant 
over the prediction horizon, 

— w M (io)~0 uj M (t) = u M (t ) V* e [t , t + T}. 

(30) 

The parameters are taken from a real machine: R = 0.86 
tt, L = 6 mH, n p = 3, K = 0.236 Vs, rated speed 314 
rad/s, rated torque 10 Nm. The trajectory optimization shall 
minimize the control error 

Pctrl(t) = (T M (t)-T* M ) 2 (31) 

and at the same time optimize the power efficiency by 



minimizing the dissipated energy 

Pioss (t) = R(i 2 d + i 2 q ) + ^ {(Li d + Kf + i 2 q ) , (32) 

where R m = 180017, thus the goal is to minimize J = 
Jo(qP c trl(t) + Pioss(t))dt + qTP ctrl (T) with the weight 
q = 20. The flat outputs id and i q are each parameterized 
with a degree 5 power series. The prediction horizon is T = 2 
ms and the sampling rate is 10 kHz. The optimization must 
also maintain the current constraints 

^ + ^</Lx = 10 2 A 2 (33) 
to protect the power inverter, and the voltage constraints 

^+« 9 2 <^a, = 330 2 V 2 (34) 

to obtain a feasible trajectory. These constraints are lin- 
earized as shown in Fig. [3] similar as in [20]. The argu- 
mentation is that only i d < makes sense (field-weakening, 
see d32l), and that the current and voltage range in q-axis is 
more important for PMSMs with isotropic rotors. 




Fig. 3. Linearized current and voltage constraints. Circle: feasible set 
of current and voltage vectors, grey: feasible set after linearization of the 
constraints. 

Numerical simulation results are shown in Fig. [4] The 
results of the QP solver are shown in blue, and the results of 
the LP solver in red. Experimental results are shown in [21]. 
At t = 0.01 s, a speed setpoint change from oj^j — to 420 
rad/s is commanded. Then, at t = 0.07 s, a load torque step 
of 8 Nm is applied. From the computational results, QP takes 
17 iterations in the worst case, wereas LP takes 24 iterations. 
This increased number is logical, as there are 20 parameters 
instead of 10. As the LP algorithm is much simpler, this 
still represents a considerable computational saving. With 
an optimized C implementation, the worst case computation 
time of the LP was 20/is on a PC (1.4 GHz clock). Runtime 
is further discussed in [21]. 

From a qualitative standpoint, the resulting behavior is 
identical for both solvers. The system is always operating 
at its performance limit. There is a (feasible) voltage peak 
on v q at t = 0.01 s to rapidly establish a torque current 
i q , and after this, the induced voltage increases proportional 
to speed on Vd and v q . The iron losses are reduced by 
imposing i d < 0. Then, at high speed, the load step 
is quickly compensated via the cascaded speed controller, 
which commands a torque increase. To do this, the trajectory 
optimization generates a negative peak on id (green arrow) 



which, according to the model, reduces the induced voltage 
on i q , thereby j^i q is higher as when keeping id small. This is 
the advantageous behavior of predictive control; in contrast 
to feedback controllers with saturation, the coupling between 
the states is exploited to bypass the constraints in an optimal 
fashion. This result is only visible as a high optimization 
horizon and the constraints are included in the trajectory 
generation. 

As the constraint handling of LP is suboptimal, there is 
a quantitative difference, especially in dynamical operation. 
The peaks are generally somewhat smaller, and at t = 0.07, 
the negative ? d -peak is of shorter duration and thus less 
effective. 

V. Conclusion 

A trajectory optimization scheme minimizing a quadratic 
cost functional to generate continuous trajectories for a 
linear control system with linear constraints was presented. 
The scheme recombines several methods in order to obtain 
a computationally efficient solution. These are the use of 
differential flatness and of a parameterization using a poly- 
nomial basis. A result on polynomial nonnegativity is used 
as a suitable way to parameterize the constraints. Further 
developments are done to formulate a linear programming 
problem. 

Numerical simulations of a predictive controller for a 
permanent magnet synchronous machine show that the sub- 
optimality of the method is acceptable. The optimization 
problem is sufficiently simplified and solvable in real-time, in 
the presented application, 2 ms prediction are implemented 
at a sampling rate of 10 kHz. 
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Appendix I 
Proof of Proposition III. 2 

Convexity of the conditioned cost function. 
Assume Q is positive definite. Then, 



x T (t)Qx(t) > 0, Vx(t) j= 



(35) 



and 



J= x T {t)Qx(t)dt > 0, x(t) 7^ OVi 6 [0, T], (36) 
Jo 

The inverse model replaces x(t) by polynomials p(t) with linear 
coefficients in a 

J = / p T (t)Qp(t)di > 0, x(t) = p(i) / (M G [0, T], (37) 



Define the primitives P(t) = / p(t)dt thus P(T) = J Q T p(t)dt 

J = P T (T)QP(T) > 0, VP(T) ^ 0. (38) 

The expression P(t) is the primitive of a polynomial, thus a 
polynomial, with linear coefficients in a, and P(T) is rewritten 
as 

P(T) = Aq (39) 
and we assume rank(A) = n with n — dim(x). It follows 

(40) 



J = a T A T QAa > 0, VP(T) ^ 0. 



The matrix of the parameterized cost function is K = A QA and 
as we assumed Q is positive definite we know Q = B T B (Cholesky 
decomposition). The weight matrix is then 



K = A T B T BA = (BA) T (BA) 



(41) 



which is positive definite as any matrix K = C T C for some C 
with rank(C) = n is positive semidefinite [CD. Meyer, Matrix 
analysis and applied linear algebra, SIAM books, 2000, pp. 566]. 
The parameterized cost functional J is thus a convex function of 
the parameters a. | 



Appendix II 
Proof of Proposition III. 3 

Sufficient affine conditions for the constraints. 
The polynomial 



P(s) = J2 ClS * ^ °> 



(42) 



with Cj G R, is analyzed on non-positivity over a segment s G [0, 1]. 
A first necessary and sufficient condition is 



P(0) = c < 



(43) 



which in the following is assumed satisfied. Furthermore, the TV 
conditions 



P UJ- 0, k=1 - N 

are also assumed to hold for all Ci. 

These conditions can be rewritten in matrix notation 

c + Qc < 



(44) 



(45) 



with c = (c ,..,c ) T G M. N , c = (c 1 ,..,c N ) T G and Q G 
R JVxAr such that 



i = 1..NJ = 1..N. 



(46) 



It can be shown that det(Q) 7^ for TV > 0, and that Q is positive 
definite. It follows 

c < -Q _1 c (47) 
which can be placed into the polynomial equation 

P(s) = co + s T c < c - s T Q _1 c = (-co)e (48) 



with s = (s, .., s N ) T € R N and e 



-1 + s T Q- 1 (l,..,l) T . As 



we assumed Co < 0, the upper bound of P(s) under the mentioned 
conditions is at when e is at its maximum. It can be shown that 
the upper bound of e, A = sup{e} Vs G [0, 1], is positive and 
only dependent on TV, as Q is known. Some values, which were 
computed numerically, are shown in the table below. 



TV 


2 


3 


4 


10 


20 


A = supje} 


0.125 


0.064 


0.041 


0.012 


0.005 



Therefore if the conditions d44b hold, we have 

P(s) < -AP(0). 



(49) 



Shifting the conditions by the constant (and negative) factor AP(0), 
the sufficient conditions for non-positivity of the polynomial P(s) 
are found: 



P(0) < 0, 
P(-|)-A-P(0) <0, k=l..N. 



(50) 
(51) 



600 




50 



-100 



0.02 



0.04 0.06 

Time t [si 



0.08 



200 



-200 



locus curve /' 
voltage /' 

I 



-500 



v d [V] 500 



800 



g ouu 
400 

T3 



200 



0.02 



0.04 0.06 

Time f [s] 



0.08 















V d 

















300 
g 200 
s 100 







0.1 





10 
5 


-5 
-10 



0.02 0.04 0.06 0.08 



0.02 



0.04 0.06 



0.08 



0.02 



0.04 0.06 

Time t \s\ 



0.08 







locus curve'' 
current / 












I 

/ 


'n 




✓ 



















p_ 










UU 


















0.1 



0.04 0.06 

Time t [s] 



0.1 



0.1 













— ^f- 






/ V 

/ q 








■ ■ ■ ■ f f 








0.1 




0.1 



Fig. 4. Results of predictive control applying the presented trajectory optimization. Red: results with LP solver. Blue: results with QP solver. 



